Generalized Gibbs ensemble and work statistics of a quenched Luttinger liquid 
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We analyze the probability distribution function (PDF) of work done on a Luttinger liquid for an 
arbitrary finite duration interaction quench and show that it can be described in terms a generalized 
Gibbs ensemble. We construct the corresponding density matrix with explicit intermode correla- 
tions, and determine the duration and interaction dependence of the probability of an adiabatic 
transition and the PDF of non-adiabatic processes. In the thermodynamic limit, the PDF of work 
exhibits a non-Gaussian maximum around the excess heat, carrying almost all spectral weight. In 
contrast, in the small system limit most spectral weight is carried by a delta peak at the energy of 
the adiabatic process, and an oscillating PDF with dips at energies commensurate to the quench 
duration and with an exponential envelope develops. Relevance to cold atom experiments is also 
discussed. 
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Introduction. Non-equilibrium many-body dynamics 
constitutes a terra incognita in comparison to its equilib- 
rium counterpart. Its exploration has begun recently by a 
series of experiments on cold atomic gases [iHj] and other 
systems Q, triggering valuable theoretical works [flQ- A 
number of interesting issues has been analyzed, such 
as thermalization and equilibration and their relation 
to integrability, defect and entropy production due to 
universal near adiabatic dynamics, quantum fluctuation 
relations[8|, non-linear response etc. 

Monitoring non-adiabatic dynamics provides a great 
deal of information about the universal features of the 
quantum system at hand. The scaling of expectation 
values or the first few moments of observables (e.g the 
defect density) after a quench through a quantum critical 
point can be expressed in terms of the equilibrium crit- 
ical exponents [6J, [7[. However, the full characterization 
of a quantum state is only possible through its all higher 
moments, encoding unique information about non-local 
correlations of arbitrary order and entanglement [9[. This 
is equivalent to determining the full distribution function 
of the quantity of interest. While its equilibrium evalua- 
tion is already rather involved [9|, obtaining the full non- 
equilibrium distribution function of a physical observable 
has rarely been carried out[10[. 

A delightful exception is the statistics of work done 
during a quench, which has been studied in Refs. [Illll2| 
for a sudden quench between gapped phases, separated by 
a quantum critical point (and gap closing). The proba- 
bility distribution function (PDF) of work done, P{W), 
involves all possible moments of energy [8[, thus providing 
us with full characterization of the energy distribution. 

While the transition between two gapped phases is of 
great interest, many interacting one-dimensional systems 
form gapless Luttinger liquid (LL) states [l3j. In partic- 
ular, interacting cold atoms in a one dimensional trap, 
e.g., often form such LL's, as also confirmed by experi- 



ments [2|, y, Il4l4l6| , but LL states appear in various spin 
models or interacting fermion systems [l3| . This state of 
matter is characterized by bosonic collective modes as el- 
ementary excitations, and by especially strong quantum 
fiuctuations. How this system reacts to a time dependent 
protocol, i.e. a quantum quench, is a highly nontrivial 
problem, though some of its properties have already been 

analyzed [il-ij- 

Here we shall study the PDF of work on this proto- 
typical example of a Luttinger liquid, determine P{W) 
after an arbitrary quench protocol, and also construct 
explicitly the generalized Gibbs ensemble which repro- 
duces all moments of P{W). We remark that this is one 
of the rare occasions, where the generalized Gibbs en- 
semble can be constructed analytically for an interacting 
model. The study of an arbitrary quench protocol is 
inspired by the observation that, in reality, quenches are 
neither completely adiabatic nor instantaneous and, — 
as we demonstrate through the properties of P{W), — 
the characteristic quench time is a crucial parameter of 
the quench itself. The work PDF is found to exhibit sev- 
eral universal forms (Gumbel or exponential distribution, 
e.g.), as controlled by the system size and interaction de- 
pendent many-body orthogonality exponent, a, and the 
duration of the quench. 

Hamiltonian. We consider an inherently gap- 

less system of hard core bosons (or an initially non- 
interacting Fermi gas) in one dimension, which is inter- 
action quenched by a given protocol into a final LL liquid 
state. The corresponding LL Hamiltonian reads 
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Here u!q{t) — v{t)\q\, and v{t) = v + Sv Q{t), with v the 
bare " sound velocity" , Sv its renormalization arising from 
interaction, and 6^ the creation operator of a bosonic 



density wave. The interaction gq{t) — g2{q)\q\Q{t) and 
the velocity are changed within a quench time r, with the 
quench protocol Q{t) satisfying Q{t < 0) = and Q{t < 
t) ~ 1. For a linear quench, in particular, Q(0 < t < t) = 
t/r. Eq. dl]) constitutes the effective model for bosons 
quenched away from the hard-core limit as well as for 
fermions quenched away from the non-interacting limit, 
or for an XXZ spin chain [13|, l20| , though our findings 
apply to interacting initial states as well [22] ■ 

Since Eq. ([T]) is quadratic, the time evolution can be 
formally determined exactly. From the Heisenberg equa- 
tion of motion, we obtain [20i 



b,it) = u,{t)b,{o) + v;it)bt,io), 



(2) 



where the time dependence is carried by the time depen- 
dent Bogoliubov coefficients Uq{t) and Vq{t), satisfying 
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with the initial condition Uq{0) = 1, Wg(0) = 0. 

Generating function of work. Armed with the for- 
mal solution of the time-dependent Bogoliubov equa- 
tions, Eq. (ID), we analyze the statistics of work done. 
Albeit the work done has been studied in classical statis- 
tical mechanics exhaustively, its quantum generalization 
has been carried out only recently [8| , and its properties 
are known for very few systems. The quantum work can- 
not be represented by a single Hermitian operator (<^ not 
an observable), but rather its characterization requires 
two successive energy measurements, one before and one 
after the time dependent protocol (thus work character- 
izes a process). The knowledge of all possible outcomes 
of such measurements yields the full probability distribu- 
tion function (PDF) of work done on the system. 

The characteristic function of work after the quench, 
G(A) = J dW e'^^ P{W) can be expressed as 

G(A, r) = {cxp[iXHH{t > r)] exp[-zAffH(0)]) , (4) 

where Hnit) is the Hamilton in the Heisenberg picture, 
and the expectation value is taken with the initial ther- 
mal state. For a sudden quench (SO), r = 0, and G(A, r) 
coincides with the Loschmidt echo (ll| . The expectation 
value, Eq. Q is independent of i for i > r, but depends 
on the quench protocol and its duration t. Hnit) is ob- 
tained by expressing the time dependent boson operators 
in Eq. ([!]) using Eq. ^. Eq. @ can then be evaluated 
at T = temperature using identities familiar from the 
theory of squeezing operators, yielding [22| 

In GiX,T) = iXEad -Y,\n{l + nq{l ^ e^'^^^)) , (5) 

9>0 

with Ead = Ef — Ei the difference between the adia- 
batic ground state energies in the final and initial state, 



and Uq = [ujq{t) - ftq + 2Im{v*{t)dtVq{t)}] /2rtq the oc- 
cupation number of mode q in the final LL state, and 

fig = \/uJq(t > t) — gq{t > t) the corresponding excita- 
tion energy [l^ . 

Generalized Gibbs ensemble. The fact that Eq. ([S]) 
depends only on the occupation numbers of the steady 
state indicates that a generalized Gibbs ensemble (GGE) 
may describe the final state [7[. The analytic construc- 
tion of the final density matrix is usually an inadmissible 
task. Therefore, one typically focuses only on few body 
observables, and tries to build an approximate density 
matrix describing these. Such an approach is, however, 
unable to account for the complete PDF of work, which 
depends on all possible moments of energy. 

In our case, the final Hamiltonian can be diago- 
nalized by a Bogoliubov transformation giving Hf = 



l^q^O "g"-? 



Ef. In the steady state (i » r), the n^'s 
and their arbitrary products are constants of motion, and 
therefore the density matrix of the GGE should be built 
up, in principle, from all of these operators [181 . We 
find, however, that for a T = temperature quench the 
density operator 
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(6) 



accounts for all intermode correlations in the final state. 
Here the mode dependent inverse temperatures /3g are 
defined through riq = {fiq) = l/[exp(/3gi7g) — 1], and 
Zg = Tr{exp[— ^^Q /3gJ7gfig]}. Indeed, it is easy to 
show that Tr{pG e*-^(^/~^^)} reproduces G(A, r) and thus 
the complete work distribution[22|. Moreover, it gives 
back the expectation value of any operators in the steady 
state. Notice that the delta-functions in pc imply perfect 
correlations between the mode pairs ±g. 

The structure of Eq. (|6]) follows from the observation 
that, while time evolution does not conserve the num- 
ber of bosons in a given pair of modes ±q, it preserves 



At) 



j{t). Since the only non-zero element of the 



initial density matrix corresponds to rig — n_g = at 
zero temperature, this can only evolve along the diagonal 



"direction", nq{t) 



j(t) = 0. The assumption that 



evolution during the quantum quench thermalizes the en- 
ergy distribution of a given momentum pair with this 
constraint then amounts in the density matrix, Eq. ([5]). 
A given pair of modes thus thermalizes only along the 



diagonal of the density matrix, fig 



-q^ 



characterized 



by an effective inverse temperature /3g, while the weight 
of the non-diagonal states rig 7^ ri_g remains zero, as in 
the initial state. Though umklapp processes may lead to 
further thermalization at larger time scales, this struc- 
ture is expected to be stable within experimental time 
scales [22 1 . 

Perturbative generating function. Though an exact 
solution is formally also possible, the general properties 
of the final work PDF are already captured by a more 



transparent pcrturbative solution of Eq. (j3]) [20[. We 
thus expand Eq. ([5]) for small 52(9) and Sv, and get for 
large system sizes L 
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Here Ead = -{L/v){g2/vTo)yi6TT + • • • < and f{t) = 
To/(t + iTo), with To a short time cut-off associated with 
the finite range of interaction, (72(9) = 52 exp(— To^lgl). 
Interestingly, the velocity rcnormalization, Sv does not 
enter to lowest order. The cumulants, C„ of the work 



done can be derived by expanding Eq. ([7]) in A (see [23 ) ■ 
Work PDF: generic properties. To analyze the PDF 
of work it is worth introducing the dimensionless work, 
measured with respect to the adiabatic ground state en- 
ergy shift. 



W=iW-Ead)/\Ead\ 



(8) 



The distribution of w is then obtained by Fourier trans- 
forming G(A, r) as 



p{w) = Vad S{w) + p{w) 



(9) 



The Dirac-delta peak corresponds to the probability of 
staying in the adiabatic ground state, while the broad 
structure p(w) is associated with transitions to excited 
states with w > 0. The weight Vad can be expressed as 

r r 
\n{Vad) = -*a / [ dtidt2Q'{h)Q'{t2)fiti-t2) . (10) 



The prefactor a = \EadTo\ ^ N{g2/v)'^ denotes the total 
angle of Bogoliubov rotations {N ~ L/vtq is the num- 
ber of particles), and can be viewed as the many-body 
orthogonality exponent. It is also closely related to the 
fidelity susceptibility |23| . Alternatively, we can rewrite 
it as a ~ L/l with I the mean free path. Thus a ^ 1 de- 
scribes, using fidelity nomenclature, the thermodynamic 
/ small system limits [23[ or, alternatively, corresponds 
to the diffusive/ballistic limits, respectively, depending 
on the picture used. 

In the adiabatic limit (r — )■ 00), a finite system al- 
ways stays in its ground state, and the time evolved wave 
function coincides with the lowest energy eigenfunction 
of the instantaneous Schrodinger equation j24| . Con- 
sequently, only the first term remains in Eq. ^ with 
Vad = 1- For T ^ To, on the other hand, Vad scales as 
^ exp(— a) ^ exp(— csi. L) (see Fig. [T|), and in the limit 
L — ^ cx) — but fixed interaction — Vad vanishes due to 
the orthogonality catastrophe. 

Sudden quench (SQ) limit. In the extreme limit 
of a SQ, T <C To, G(A, t) simplifies to G(A) ~ 



exp [i_BadA^/(A -I- iTo)] , and the continuum part of the 
PDF of work is evaluated exactly as 

Psq{w) = Vad exp(-a'u;) a w'^l"^ h{2ayfw) , (11) 

with Vad = exp(— a) and I\(x) the modified Bessel func- 
tion of the first kind. This is the non-central x^ distribu- 
tion with non-centrality parameter 4a in the limit of zero 
degrees of freedom [25|. The average work is zero [20|, 
since for a SQ the system remains in its initial state and 
— on average — there is no back reaction. Entropy is, 
however, generated by populating high and low energy 
configurations. 

The shape of piw) depends crucially on the orthogo- 
nality parameter, a. In the thermodynamic limit, a 3> 1, 
almost all probability weight is carried by a peak centered 
at around W = ^{w = \) and of width ^W ~ \Ead\l\fo., 



PSQ («^»""') 



exp I —a [1 — a/w] 



(;3/4v^ 



(12) 



whose high energy tail decays according to the Gamma 
distribution, ^ exp(— aw)/'u;'^''*. In the small system 
regime a <C 1, on the other hand, the delta function re- 
tains almost all weight, and transfers only a fraction ^ a 
to an exponential distribution of width AVF ^ \Ead\loL 
and threshold at Ead for w <C ot"^ . In the cross-over 
regime, a ^ 1, the maximum shifts to lower energies and 
the PDF of work develops a sizable value right above the 
threshold at Ead (see Fig. [T]). The maximum of P(W) 
occurs zXW> Ead for a > 2, while the PDF becomes 
monotonically decreasing for a < 2. 

Finite quench times. For finite duration quenches, 
in addition to the orthogonality parameter a, the work 
statistics also depends on t and the protocol Q{t) itself. 
For dcfiniteness, we focus here on a linear quench |22l |. 
and measure the degree of adiabaticity by f = t/to . 

For a finite duration quench, t > 1, only a fraction 
1/f of the excitations experiences the quench as sudden. 
Consequently, in the expression of T^ad, the orthogonality 
exponent a is replaced by Ur ^ ct/r, and Vad becomes a 
monotonously increasing function of t (see Fig. [T]) . The 
crossover with increasing ctr from Vad ^ 1 to vanishingly 
small spectral weight, Vad, occurs at a '^ t. 

Close to the threshold, W — Ead "C 1/t, only states 
with energy smaller than 1/t and thus feeling a SQ con- 
tribute to work. Therefore, apart from a normalization 
factor, the PDF of work agrees with the SQ result. 



p(w < a ^f ^) fiiVadexp{a)psQiw) 



(13) 



and depends on t only through Vad- 

For f ^ 1, however, Eq. (fT3)) describes only a small 
region close to Ead (see thin black lines in Fig. [TJ, 
and the overall shape depends both on a and f. For 
4a ^ f, almost all spectral weight is carried by the 
non-adiabatic processes (p{w)) around the typical value 



a = 20 (thermodynamic limit) 



a = 4 (crossover region) 



a = 0.2 (small system limit) 

0.04 /lA 




-0.5 



FIG. 1. (Color online) The PDF of work done on a LL is plotted after a linear quench from the numerical evaluation of Eq. 
(blue solid line). Left panel: a = 20 with f = 0, 1, 2.5 and 5 from right to left and 180 (inset, P{W > Ead) only); middle 
panel: a = 4 with f = 0, 1, 2 and 4 with increasing peak height and 55 (inset); right panel: a — 0.2 with f = 0, 2, 5, and 25 
from right to left. The thick magenta line denotes the exact SQ expression (Eq. (|lip ). the red dashed line represent Eq. (|15|l . 
the thin black line in the middle panel visualizes Eq. (|13p . while the green dash-dotted line shows Eq. (|14|l . The vertical arrow 
a.t W = Ead denotes the Dirac-delta peak, whose spectral weight Vad is shown in the inset of the right panel on semilog scale 
as a function of the ramp time r. 



Wtyp — Ead "^ 2|£'ad| ln(f)/f^, clearly separated from 
the adiabatic process. For f ^ 4a, the adiabatic pro- 
cess gains spectral weight, Vad ~ 1, but a maximum for 
W > Ead remains present, though it gradually merges 
with the adiabatic processes. 

In particular, in the small system limit Ot- <C 1, e.g., 
we can expand Eq. ([7]) to get 



2 fsm\waf/2\ 

p{w) « Vada z-r— 

' waT/2 



exp(— aui) . (14) 



Thus, for a linear quench, for energies commensurate 
with the quench time the PDF is zero (see Fig. [1]). This 
is related to the steady state behavior of the occupation 
numbers, n^ « [(72(9) sin('i;|(7|r)/2w^|(7|r]^ for g2 ^ v j22l |. 
reflecting that modes with energy commensurate to the 
quench time stay almost unoccupied at T = 0. In 
this limit (a <C r and f ^ 1), the system evolves almost 
adiabatically, non-adiabatic processes have only a small 
probability ^ a/f, and the typical work done in case of 
a rare non-adiabatic process is Wtyp ~ —o^tt/t. 

Increasing a, the zeros of the PDF turn gradually into 
dips, and the PDF develops a more universal form. In 
the thermodynamic limit ar ^ I, using the method of 
steepest descent we obtain 



p(w) K , I cxp w , 

^ ' 2Vtan3(s)7r V V 2 



(15) 



for as ^ f, with s = arctan[-\/exp(wf2) — 1]. For 
w :» l/f^ >• l/a^, p(w) in Eq. (fTS)) behaves as a gener- 
alized Gumbel distribution of index a = ^ -I- |f J26J . 
This latter emerges in the context of global fluctua- 
tions, describing the limit distribution of the a-th max- 
imum of a sequence of independent and identically dis- 
tributed random variables 9|. The distribution in the 



l/f2 ^ w ^ Xjo? region resembles closely to Eq. 
apart from its normalization. 

Experimental relevance. Our results can be tested on 
one-dimensional hard-core bosons |27| or non-interacting 
fcrmions as initial states. The detection of the PDF 
of work requires two energy measurements, one before 
and one after the time dependent protocol. The first 
energy measurement can be omitted if we prepare the 
initial wave function in an energy eigenstate of Hit = 0). 
The resulting energy distribution can then be probed 
using timc-of-flight experiments [j, |7| , similarly to Ref. 
28|. The crossover between the various regimes can be 



monitored by tuning t/tq and a ^ N ((72/^) , where 
N is the number of particles in a ID trap, typically 
with iV - 10^ - 10^ atoms 0, S, [H. By choosing 
92/ V ~ l/viV, a becomes of order unity, facilitating the 
observation of crossover between the various regimes. For 
one-dimensional interacting bosons (i.e. Bose-Hubbard 
model), V ^ J and 92 ^ J^ /U ior U ^ J (close to the 
hard-core boson limit) with U the on-site interaction [29| 
and J the hopping amplitude. By quenching away from 
the initial t/ ^ J <^ 32 ~ limit (e.g. by changing the 
lattice parameters or tuning the Feshbach resonance), a 
final interaction U ~ J\/N is reachable. For weakly in- 
teracting fermions, v ^ J and 92 ^ U , therefore ramping 
from the weakly interacting case to U ^ J/vN is desir- 
able. Nonetheless, our results apply also to interacting 



initial states 22 



Summary. We have studied the PDF of work done on 
a LL after an interaction quench, realizable in strongly 
interacting Bose systems. We have constructed the den- 
sity matrix of the generalized Gibbs ensemble with in- 
termode correlations, describing arbitrary correlations of 
the steady state, thus the PDF of work. The PDF ex- 
hibits markedly different characteristics depending on the 



system size, quench duration and interaction strength. 
Our method is apphcable to the full PDF of other ob- 
servables as well, e.g. density fluctuations [30] ■ We also 
emphasize that our results in Eqs. (O and 1^ apply also 
to a variety of other systems with effective bosonic Hamil- 
tonians as in Eq. ([T|) , including interacting higher dimen- 
sional bosons or spin systems within a spin- wave theory. 
This research has been supported by the Hungar- 
ian Scientific Research Funds Nos. K72613, K73361, 
K101244, CNK80991, TAMOP-4.2.1/B-09/1/KMR- 
2010-0002 and by the Bolyai program of the Hungarian 
Academy of Sciences. 
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SUPPLEMENTARY MATERIAL FOR 

"GENERALIZED GIBBS ENSEMBLE AND 

WORK STATISTICS OF A QUENCHED 

LUTTINGER LIQUID" 

QUENCHING BETWEEN INTERACTING 
LUTTINGER LIQUIDS 

We demonstrate here that our results for the PDF of 
work applies also when we quench between interacting 
initial and final states. More precisely, we do not need 
to start directly from the hard core boson limit for a 
bosonic LL or from strictly non-interacting fermions for 
a fermionic LL. Let's consider an initially interacting LL, 
given by 



9a 



H,^J2^,b+b, + '-^[b,b 



bp^,] 



(SI) 



q^O 



where g* is the initial interaction. This can conveniently 
be diagonalized by a standard, time independent Bogoli- 
ubov transformation as 



bg = cosh{(j>g)ag - sinh(0,)a+q, 
tanh(2(/),) = ^, 



yielding 






UJa 



Ulqa+ttq, 



(S2) 

(S3) 
(S4) 



where ujg = J{ijjqY ~ (ffg)^i 3'iid Ei is the ground state 
energy of Hi with respect to the non-interacting ground 
states. The ground state of this Hamiltonian is the vac- 
uum of the a bosons. 

The interaction quench is described, in the language of 
the initial bosonic operators, by the additional term 



H' 



Agg{t) 



9#0 



bqb^q 



Kb^-,] 



(S5) 



where Ag^it) = {g^ ~ gV)Q{t), and g' is the final interac- 
tion strength. After applying Eq. (|S2|) to this, we obtain 
the total time dependent Hamiltonian, H = Hi + H' as 




^9qit) 9\ 



OJn 



a^aq+ 



+ ^^Ka_, + «>!,] 



^9qit) 9l 



(S6) 



initial and final interactions, provided that we stay in 
the perturbative regime throughout the quench, namely 
ujq ^ ISqU^II- The threshold above which work is pos- 
sible, occurs at the difference of the adiabatic ground 
state energies of the final and initial states, namely at 
Ead = -{L/v)[{gif - {g\f]/{vT^fl&n +.... The or- 
thogonality exponent is determined by another energy 
scale, Eoe = \Ead{gi - gh)/{g2 + 52)l> giving a = EoeTo- 
After redefining w = {W — Ead)/ Eoe, our results for p(w) 
hold for this general case as well. 



CALCULATING THE CHARACTERISTIC 
FUNCTION OF WORK 

The Hamiltonian in the Heisenberg picture for i > r 
is given by 

HH{t) = Y, Co{q, t) {b+bq + b^qbtq) " LOq{t) + 
q>0 



+ci (<7, t)bqb-q + c* {q, t)b+b 



+ 

q "-q- 



(S7) 



where 



Coiq,t)^^q{t){\Uqit)\^ + \Vq{t)\^) + 

+gq{t)2Rc[uqit)v;it)], (S8) 

Cliq,t)=2uJq{t)Uqit)Vq{t)+gqit)iul{t)+vlit)). (S9) 



This is simplified upon realizing that the operators 



Ko{q) 
K+{q)^b+b+_q, 



b+bq + b^qbtq 



K^{q)^bqb.q 



(SIO) 
(Sll) 



are the generators of a SU(1,1) Lie algebra, satisfying 
[K+{q),K_[q)] = -2Ko[q), [Ko{q), K±{q)] = ±K±{q), 
and the operators for distinct g's commute with each 
other. Following Ref. [31|, we obtain 

cxp[iXHH(t)] = Y[cxp[(3+{q,X,t)K+{q)-iLjq{t)X\x 
q>a 
X exp[2/3o(g. A, t)Ko{q)] exp[/3_ (q. A, t)K_{q)l (S12) 

and the /3±^o{q, A, t) coefficients can be determined using 
Ref. 3l|. When starting from the ground state at T = 0, 
the characteristic function of work is obtained as 



G(A, r) = exp [ ^ Mq, A, t) ~ iUq{t)\ j 

\9>0 / 



(S13) 



with 



This, apart from the constant first and last terms on the 
r.h.s, is identical to Eq. (1) in the main text, after re- 
defining its parameters, w, 8v and gq{t) appropriately. 
Therefore, our results for the PDF of work apply for any 



/3o(g, A, i) = - In (cos(rjqA) 
,ujq{t) + 2lm[vl{t)dtVq{t)] 



n„ 



smiriqX) , 



(S14) 



where ftq = , /w^(i > r) — g'^{t > t) is the time indepen- 
dent adiabatic eigenencrgy of a given mode in the final LL 



state 
as 



13| . After some trivial steps, Eq. (|S13p is rewritten 



In (g{X, t)] =-Y^ In (1 + 71,(1 - exp[2ingX])) 



q>0 



~i\E, 



adi 



where 



UJq{t) +2lm[vl{t)dtVq{t)] 
2Q.a 



is the occupation number in the steady state and 

q>0 



(S15) 



(S16) 



(S17) 



is time independent, and stands for the difference of the 
adiabatic ground state energies of the final and initial 
states. 

The occupation number in the steady state is obtained 
for a linear quench from Eq. (|S16[) using the results of 
Ref. in as 



52(g) sin(w|g|r) 



2^2 



191^ 



(S18) 



for g2 ^ V. We have checked this prediction by numer- 
ically solving the differential equation in Eq. (3) in the 
main text. The occupation numbers turn out to be peri- 
odic in the time averaged mode energies, defined by 



(q) = - / ^oj^q{t)-g^q{t)dt, 



(S19) 



and the numerically evaluated occupation numbers are 
plotted in Fig. ISII States with average energy commen- 
surate with T possess the lowest effective temperatures. 
The zeros predicted by Eq. (JSISP turn to sharp dips with 
increasing interaction, and our analytical expression de- 
scribes rather reliably the numerical data. Although for 
finite 32, all effective temperatures are finite, the steady 
state is far from being thermal, as evidenced by the highly 
non-thermal structure of the steady state density matrix. 



ON THE STEADY STATE DENSITY MATRIX 
FROM GGE 

The steady state density matrix reveals a highly non- 
thermal structure, namely 



PG = Y~,'n. '^^P [-^<3^<3"-<?] ^n,,n- 



(S20) 




12 3 4 5 

iOav{q)r /2tv 

FIG. SI. The numerically evaluated occupation numbers in 
the steady state are plotted on a semilogarithmic scale after 
a linear quench with 5v — 0, (72 /« = 0.1 (blue solid line), 0.5 
(red dashed line) and 0.8 (black dash-dotted line), ranging 
from weak to strong interactions. The result of Eq. (|S18|) 
is practically indistinguishable from the blue solid curve ex- 
cept close to ujav{q) = 2im/T with n integer, where the ef- 
fective temperatures are the lowest. The figure is valid for 
arbitrary q since it depends only on the dimensionless combi- 
nation LiJav{q)T. 



possessing finite matrix elements only along fiq — h-q. 
Thermalization would mean the same matrix elements for 
a fixed hq -f fi_g, which is not the case here. Therefore, 
one can ask to what extent this density matrix is immune 
to additional perturbations, not considered within our 
Hamiltonian. 

On the one hand, one can argue that experimen- 
tal results on one dimensional cold atoms did not find 
any sign of thermalization, but rather prethermalization, 
i.e. reaching a certain non-thermal steady state, took 
place @,y,llfl,ll5l- These experimental results are nicely 
accounted for by a simple Gaussian model, similar to our 
Eq. (1) in the main text, without additional terms. On 
the other hand, from a theoretical point of view, addi- 
tional interactions between the h bosons are generated by 
the non-linearity of the non-interacting dispersion rela- 
tion. As was investigated in Ref. 



32 



g>0 



, this broadens 
the otherwise sharp peak around the bare dispersion at 
Ljq = v\q\ in the spectral function of the bosons, which, 
as a rough estimate, scales with ~ q^ jm^ where m is 
the effective mass arising from curvature effects. In the 
long wavelength limit (g ~ 0), this gives a "lifetime", 
which is usually much longer that the typical experimen- 
tal timescales. Higher order terms in the h bosons also 
arise from density-density interactions [29J, but these pro- 
vide higher powers of q in the broadening of the bosonic 
mode, suppressing further their effect. 

Nonlinear terms of sine-Gordon type [13[ are ab- 
sent without any lattice, i.e. for interacting parti- 
cles in the continuum limit, which can also be realized 
experimentally J27| . In the presence of an optical lattice, 
these are inevitably present, though their effect can be 
weakened by choosing incommensurate fillings or sup- 
pressing spin backscattcring (e.g. by using single com- 



poncnt bosons). Therefore, our non-thermal density ma- 
trix in Eq. (|S20[) is expected to describe fairly reliably 
the steady state of interaction quenched one dimensional 
systems. 

Having established the validity of our density matrix, 
we now turn to the derivation of results, presented in the 
main text, using the steady state density matrix. The 
partition function of the GGE is determined as 



Zg = 



n i 

J-J- l-expf 



g>0 



- exp[-/?qrig] 



(S21) 



and the characteristic function of work reads as 
G'(A,T)exp(-iA£'ad) = 

_. oo 

= — J]^ ^ exp \-^qVtq -I- 2iAf2g?7,] = 
^ g>On,=0 



exp - ^ In 



q>0 



1 - exp[-/3gr^g] 
- cxp[— /3qJ7g -|- 2iAf2g] 

/exp[/?qSlq] — exp[2iAr2g 
V cxp[/3,17q] - 1 



V ^ V exp(/3,i7,)-i;y 



THE CUMULANTS OF ENERGY 

The cumulants, C„ of the PDF of work done can be de- 
rived after expanding the characteristic function of work 
done in power series as In (G'(A,t) j = Y^=\ Cn(*A)"'/n!, 
yielding 



Ead 



= 5n,l - 



^■('■)0'('^)-g,.M.. (S23, 



[i(tl - t2) - To] 



where the first term denotes the adiabatic ground state 
energy difference between the initial and final state, while 
the second one stems from the non-adiabatic evolution 
(i.e. heat). The first cumulant, Ci = {H{t > r)) was 
already calculated in Ref. |20l Their behaviour is illus- 
trated in Fig. [S2] after a linear quench: Ci — Ead decays 
as \n(TjTo)/T'^, since Q'{t) exhibits kinks at t = and 
r[20|, |34|, consequently all C„>i decay as t~^. 

For a SQ, the cumulants are Ci = {H{t > 0)) = 
EsQi= to second order in (72 and Sv within our scheme), 
C„>i = an\/TQ~'^. 
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FIG. S2. (Color online) Several cumulants of the work done 
on a LL are log-log plotted as a function of the quench time 
for a linear protocol. Close to the SQ limit (r <C tq), all 
properly normalized cumulants are equal, while in the near 
adiabatic limit (r ^ ro), these approach 2/n(n — 1). 



ANALYTICAL RESULTS FOR A LINEAR 
QUENCH 

In the case of a linear quench, the characteristic func- 
tion of work is obtained as 



In 



G(A,t) 



iEn 



A - 4 {g{\ + zro) - .g(iro)) 



(S24) 



where (7 (z) = (z + t) ln(z + T) + (z — t) ln(z — t) — 2zln(z). 
The cumulants are obtained from a variant of Eq. (IS23[) 



as 



Cn . , ,,„ 2 9"-^ A£; 



-rr- == Ki - (-1) To 

Hind 



di 



EadTo 



where the heating, AE is 



AE 



-E„ 



In 1 



(S25) 



(S26) 



and Eq. (|S25p holds true for an arbitrary quench proto- 
col. The asymptotic behaviour of Eq. (|S24[) for A — >■ 00 
gives 



G{X — > 00, t) ~ exp ( iXEad — 2— arctan (7 

X (1 + fY'^ I "^P {'^^'^^ -t) ^ ^ °°' (S27) 
[ exp (iXEad - a) r ^ 0. 

From this, the asymptotic behaviour of Pad follows as 



Pad 



exp —a 



1 



6 



(S28) 



